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I.  INTRODUCTION 


In  a  recent  report  [1]  WKB  formalism  developed  by  Budden  [2]  for  an 
anisotropic  medium  was  used  to  speed  up  calculation  of  ELF/VLF  modal  height 
gains  in  the  ionosphere  to  heights  well  in  excess  of  the  peak  of  the  F  layer. 
Justification  for  use  of  the  WKB  formalism  lies  in  the  fact  that  for  altitudes 
sufficiently  high  in  the  ionosphere  (i.e.,  >  110  km  for  VLF,  >  150  km  for  ELF 
under  daytime  ionospheres  and  >  250  km  for  ELF  under  nighttime  ionospheres) 
the  electromagnetic  field  is  in  the  outgoing  whistler  mode  and  the  ionospheric 
gradients  are  sufficiently  weak.  In  this  report  the  WKB  method  is  appended  to 
an  earlier  program  [3]  for  calculating  excitation  of  the  earth- ionosphere 
waveguide  by  point  dipoles  at  satellite  altitudes.  As  with  the  latter 
program,  provision  is  made  for  calculating  excitation  factors  for  all  electric 
field  components  generated  by  both  electric  and  magnetic  point  dipoles  of 
arbitrary  orientation.  It  can  be  expected  that  the  present  program  will 
reduce  the  computer  cost  (Uni vac  1100)  by  about  a  factor  of  1/3  in  the  lower 
ELF  band  (^75  Hz)  and  by  more  than  an  order  of  magnitude  in  the  VLF  band  when 
the  source  altitudes  are  >  500  km.  We  summarize  in  the  following  paragraphs 
the  features  necessary  for  an  understanding  of  the  present  work. 

All  height  dependent  input  quantities  are  referenced  to  a  Cartesian 
coordinate  system  (x,  y,  z)  such  that  the  z  axis  is  taken  along  the  vertical 
and  positive  into  the  ionosphere.  The  (x,  y,  z)  coordinate  system  is  referred 
to  in  this  report  as  the  input  coordinate  system.  The  transmitter  and 
receiver  are  located  in  the  x-z  plane  with  z  =  0  being  the  ground.  The  source 
is  taken  to  be  imbedded  in  a  semi-infinite,  homogeneous,  but  anisotropic  slab 
with  lower  boundary  termed  WKBHT  in  the  input  coordinate  system.  Thus, 
reflections  from  above  the  source  are  ignored.  The  primary  plane  wave  field 
radiated  by  the  elevated  antenna  at  an  earth  ionosphere  eigenangle,  0n>  is 
decomposed  into  downcoming  magneto- ionic  components  and  evaluated  at  z  =  WKBHT 
in  the  input  coordinate  system.  Next,  two  independent  full  wave  starting 
solutions  consistent  with  the  condition  of  downcoming  waves  are  assumed  at  an 
altitude  termed  LWSTHT  in  the  input  coordinate  system.  LWSTHT  must  be  in  the 
isotropic  space  between  the  ground  and  the  lowest  level  of  the  ionosphere. 
The  two  independent  solutions  are  numerically  integrated  upward  to  a  level 
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termed  TOPHT  in  the  input  coordinate  system.  At  TOPHT  the  two  solutions  are 
combined  to  yield  the  proper  modal  polarization  and  decomposed  into  the  down¬ 
coming  whistler  component.  The  full  wave  whistler  component  is  then  matched 
to  the  corresponding  WKB  component  and  the  solution  propagated  to  WKBHT  via 
the  WKB  formulas.  Finally,  the  strength  of  the  full  wave  solution  is  deter¬ 
mined  by  matching  to  the  source  generated  whistler  component  previously 
determined  at  WKBHT. 

The  present  program  is  a  variant  of  a  fields  program  [4]  which 
necessitated  starting  the  full  wave  integrations  with  outgoing  wave  and 
subsequent  downward  integration.  That  is  in  contrast  with  the  present 
requirement  of  starting  the  full  wave  integrations  with  down-going  wave  and 
subsequent  upward  integration.  To  achieve  this  modification  the  specie 
profile  and  collision  frequency  read  in  (as  was  done  in  reference  [3])  is  left 
unaltered  and  coordinate  transformations  effected  internally  in  the  program. 
That  is,  read  in  of  the  latter  quantities  is  referenced  to  what  has  been 
referred  to  previously  as  the  input  coordinate  system  (z  axis  positive  into 
ionosphere,  z  =  0  at  the  ground).  The  profiles  are  read  in  from  WKBHT  down  to 
level  LWSTHT  (or  equivalently  D)  which  is  often  taken  to  be  ground  (LWSTHT 
=  0).  The  coordinate  system  and  profiles  are  then  transformed  by  rotating 
180°  about  the  x-axis  and  translating  in  the  z  direction  by  the  amount  WKBHT. 
This  is  referred  to  as  the  prime  coordinate  system  and  represented  by  (x1  ,  y1  , 
z '  ) .  Mathematical  ly  the  transformation  is  z‘  =  -z  +  WKBHT.  This  has  the 
effect  that  TOPHT1  =  -LWSTHT  +  WKBHT,  LWSTHT1  =  -TOPHT  +  WKBHT  and  WKBHT1  =  0, 
where  the  left  hand  side  of  the  equations  are  (apart  from  the  1  )  the  program 
names  after  the  coordinate  transformation.  The  transformation  from  the  (x,  y, 
z)  system  to  the  (x1 ,  y1 ,  z1)  system  is  effected  internally  in  the  program  and 
achieves  the  goal  that  the  starting  full-wave  solutions  in  the  transformed 
system  correspond  to  upgoing  wave  and  that  the  full  wave  integrations  proceed 
as  in  the  original  full-wave  program  [4],  Transformation  to  the  primed  system 
is  achieved  in  the  subroutine  XINPUT. 

Basic  inputs  to  the  present  program  must  be  supplied  from  an  ELF/VLF 
earth- i onosphere  waveguide  program  such  as  that  of  reference  [5],  The  inputs 
include  mode  eigenangles,  derivatives  of  the  modal  function  evaluated  at  the 
eigenangles  (or  related  quantities)  and  the  elements  of  the  plane  wave 
ionospheric  reflection  matrix  as  defined  in  reference  [2]. 


For  completeness,  Budden's  WKB  formulas  [2]  are  summarized  in  section  II. 
Formulas  used  in  matching  the  full  wave  and  WKB  field  are  given  in  section 
III.  Principal  outputs  of  the  program  are  the  excitation  factors  defined  in 
section  IV.  Flow  of  the  program  and  subroutine  descriptions  are  given  in 
section  V.  Section  VI  contains  a  discussion  of  the  card  deck  arrangement  for 
input  along  with  discussion  of  a  sample  output.  Application  of  the  excitation 
factors  for  calculating  mode  sums  is  discussed  in  section  VII  where  sample 
mode  sum  plots  are  included.  Appendix  A  contains  a  listing  of  the  excitation 
factor  program.  Appendix  B  contains  a  listing  of  the  mode  sum  program  used  to 
generate  the  results  shown  in  section  VII. 


II.  SUMMARY  OF  WKB  FORMULAS 

For  plane  wave  incidence  of  an  rf  wave  on  the  ionosphere,  Maxwell's 
equations  can  be  written  as  [6] 

e1  =  -iTe  (1) 

where  the  prime  denotes  (k  ^8/3z),  with  k  the  free  space  wave  number,  T  a  4  x 
4  matrix  given  by  Budden  [2],  and  e  a  column  matrix  composed  of  components  of 
the  electric  (£)  and  magnetic  (H)  fields  of  the  rf  wave.  The  transpose  of  e 


e'  =  (E  ,  -  E  ,  Z  H  ,  Z  H  )  (2) 

x  y  ox  oy 

where  Zq  is  the  free  space  impedance.  Flenceforth  the  notation  = 

Z  fiy  and  H  =  Z  H  will  be  used. 
o  J  z  0  1 

The  matrix  T  has  four  characteristic  roots  or  eigenvalues,  (i  = 
1,2, 3, 4),  which  satisfy  the  characteristic  equation 


det  (T  -  ql)  =  0 


where  I  is  the  unit  4x4  matrix.  This  equation  is  the  Booker  quartic. 
Corresponding  to  any  root  q.  there  is  an  eigencolumn  P  =  s^1^  which  satisfies 

T  S(l)  =  qis(l)  (4) 


Let  S  be  the  4x4  matrix  whose  columns  are  sv 
nonsingular,  the  column  matrix  f  can  be  defined  as 


e=Sforf=Se 


and  it  can  be  shown  [2]  that  the  elements  of  f  satisfy 


For  points  where  S  is 


f,'  =  *iq.fi.  +  21.  .f. 


;j,k  =  1,4 


where 


r  =  -s_1s' 


The  preceding  transformation  can  be  carried  out  only  when  S  is 
nonsingular.  When  two  roots  of  the  Booker  quartic  are  equal,  two  of  the 
columns  are  usually  multiples  of  each  other,  and  then  S  is  singular. 
Near  such  points  some  of  the  coupling  coefficients  r^.  are  very  large  and  the 
points  may  be  points  of  reflection  or  points  where  coupling  between  two  up- 
going  (or  downgoing)  waves  is  very  strong.  The  present  program  cannot  be  used 
in  such  circumstances. 

When  the  species  densities  and  collision  frequencies  vary  slowly  enough 
with  height  and  where  no  two  of  the  q^  become  nearly  equal,  the  terms  of  F  are 
small  quantities  and  there  is  an  approximate  solution  for  which  the  non¬ 
diagonal  elements  of  T  are  ignored.  This  solution  is  associated  with  one 
particular  root  q.  of  the  Booker  quartic.  It  is  [2] 

J 

f j  =  exp  (“i kjzq^.dz  +  kfzr„dz)  (8) 
and  the  corresponding  field  components  [in  Budden's  notation]  are 


(E  ,  E  ,  H  ,  H  )  =  (A  .F . )  1/z  (a^q.  +  a-,  -A.,  q.A.,  a,q.  +  ac) 
x>  y’  x’  y'  j  y  3Mj  4’  j’  j’  5Hj  6' 

x  exp  (-ik/zqj.dz  +  kfVjjdz) 


where 


al  =  '(T11  +  V 


a2  T11  T44  "  T14  T41 


3 4  T14  T42  "  T12  T44 


a5  T42 


a6  T41  T12  '  T11  T42 


=  qj  +  alqj  +  a2 

=  2q,A.  ♦  (q2.  -  T.„)  (2q,  +  a, )  -  (a,b,  +  b,a,) 


2rjj  =  <AjFj>  *iqj(a3b5  '  aib5  *  a5b3  '  a5b3> 


*  a4b6  ’  a4b6  +  a6b4  '  a6b4  +  qj^a3b6  ~  a3b6  +  a4b5  '  a4b5  +  a5b4 


-  a‘b4  +  agb^  -  a^b3) 


In  these  equations  the  T..'s  are  the  elements  of  the  T  matrix  given  by  Budden 
p  389. 

I  he  essence  of  the  program  documented  in  this  report  is  the  implementa¬ 
tion  of  equation  (9)  for  altitudes  exceeding  a  height  termed  TOPHT  in  the 
input  coordinate  system.  The  mode  extracted  from  the  magneto-ionic  set  is  the 
least  attenuated  downcoming  wave.  Runge-Kutta  integrations  are  used  to 
calculate  two  independent  full  wave  solutions  at  TOPHT  from  starting  condi¬ 
tions  at  LWSTHT  in  the  input  coordinate  system.  The  full-wave  solutions  at 
TOPHT  are  combined  to  satisfy  the  mode  polarization  condition  and  decomposed 
into  magneto  ionic  components.  The  downcoming  whistler  component  is  then 
matched  to  the  corresponding  WKB  component  and  the  fields  carried  to  WKBHT  via 
equation  (9). 


III.  FIELD  MATCHING 


As  described  in  the  introduction,  it  is  necessary  in  the  primed 
coordinate  system  to  integrate  two  independent  field  solutions  from  the  level 
z'  =  TOPHT*  =  -LWSTHT  +  WKBHT  to  the  level  z‘ =  LWSTHT 1  =  -TOPHT  +  WKBHT.  The 
initial  solution  vectors  correspond  in  the  prime  system  (at  level  TOPHT1)  to 
outgoing  waves  in  vacuum  and  are  taken  to  correspond  to  vertical  and  hori¬ 
zontal  polarizations.  The  initial  solution  vectors  are  calculated  in  the 
subroutine  WF  INIT  and  are  denoted  by  the  matrix  P  as  shown  below 


P(l.l)  =  Ex  =  cos(6n)  =  C. 


P(l,2)  =  Ex  =  0. 


P(2,l)  =  -E  =  0. 


P(3,l)  =  Hx  =  0. 


P(2,2)  =  -E  =  1. 


P(3,2)  =  Hx  =  cos(0n)  =  C. 


(15) 


P(4 , 1)  =  H  =  1. 


P(4,2)  =  Hy  =  0. 


where  E.  and  H.  represent  the  jth  component  of  the  electric  and  magnetic 

J  J 

fields  of  the  rf  wave  in  the  prime  system  and  6  represents  the  eigenangle  for 
t  h 

the  n  earth  ionosphere  waveguide  mode.  The  initial  solutions  are  integrated 
to  the  level  z1  =  LWSTHT'  via  the  Runge-Kutta  integration  scheme  implemented 
in  the  subroutine  WF  INTG.  At  z'  =  LWSTHT1  each  independent  full-wave 
solution,  P.  is  decomposed  into  four  magneto-ionic  components  as  follows 


P.  =  I  P(i,j)  S ( j ) ;  i  =  1,2 
j  =  1 

where  S(j)  are  the  magneto-ionic  solution  vectors.  Thus, 


(16) 


PC i  ,j)  =  S  1{ j)P. 


(17) 


where 


S_1(j)S(k)  =  0  j*k 
=  1  j^k 


(18) 
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Because  of  the  sorting  achieved  in  the  subroutine  WF  SORT,  the  upgoing  (in 
the  prime  system)  whistler  mode  is  S(2).  The  total  whistler  component  is 
structured  from  a  linear  combination  of  the  components  associated  with  the  two 
independent  full  wave  solutions.  The  two  independent  solutions  are  combined 
in  such  a  way  that  they  satisfy  the  modal  polarization  condition  at  LWSTHT  in 
the  input  coordinate  system 


Ey  1  "  iiRi,iiRh 

y  j.RiuRi 


1  -  j.Ri.i.R.1 


In  this  equation  R  is  the  plane  wave  reflection  coefficient  associated  with 
everything  above  LWSTHT  and  R  is  the  plane  wave  reflection  coefficient  from 
everything  below  LWSTHT  in  the  input  coordinate  system.  Consistent  with  con¬ 
vention,  the  first  subscript  on  R  or  R  applies  to  the  polarization  of  the 
incident  wave  whereas  the  second  subscript  applies  to  the  polarization  of  the 
reflected  wave,  with  ,,  denoting  vertical  polarization  and  x  denoting  hori¬ 
zontal  polarization. 

The  total  upgoing  (in  the  prime  system)  whistler  component,  Sy(2), 
obtained  from  decompos i tion  of  the  two  full  wave  solutions  at  z1  =  LWSTHT',  is 
then  given  by 

ST(2)  =  ((S<1,2)  -  (3(2,2) f)S(2)  (20) 

Note  that  the  minus  sign  before  (3(2,2)  follows  from  the  fact  that  P(2,2)  =  -E 
(see  equation  (15)).  Equation  (9)  with  j  =  2  is  normalized  such  that  it 

equals  equation  (20)  at  z1  =  LWSTHT'  and  the  solution  propagated  to  z'  = 
WKBHT 1  =  0  via  the  WKB  formula.  If  the  solution  at  z'  =  0  is  denoted  by 


0  via  the  WKB  formula.  If  the  solution  at  z 


[((3(1,2)  -  (3(2,2)f)S(2)]  ,_g,  then  the  strength  and  final  unknown,  A^,  of  the 
full  wave  solution  is  determined  by  [7,  8] 


BP  (1 ,2)elkQ(2)Zo 


a  -  _ 111  _  (Pi) 

m  "  (Q(2)  -  Q(1))(Q(2)  -  Q(3))(Q(2)  -  Q(  4)  [  ( (3(  1 , 2 )  -  (3(2 ,2)f  )S(  1 ,2)  ]z ,  =Q  ^ 

where  S( I , 2 ) ( I  =  1,2, 3, 4)  represents  the  1^  field  component  of  the  magneto¬ 
ionic  mode  (i.e.,  I  =  I  *  E.  I  =  2  -  -E  ,,  I  -  3  *  H  .  1  =  4  •>  H  ).  Any  of 


the  field  components  can  be  used  to  determine  A^.  The  program  uses  or  I 
=  4.  The  free  space  wave  number  is  denoted  by  k,  the  Q's  are  Booker  quartic 
solutions,  and  is  the  transmitter  location  in  the  prime  coordinate  system. 
BPm(I,J)  represents  an  amplitude  factor  associated  with  the  decomposition  of 
the  primary  source  into  an  upgoing  (in  the  prime  system)  magneto-ionic 
component  (J  =  2)  associated  with  the  eigenvalue  0  .  It  is  determined  from 
equations  (17),  (22),  (24),  (25),  and  (56)  of  reference  [8]  and  will  not  be 
reproduced  here.  The  subscript  m  on  A  and  BP  denotes  the  source  type  and 
configuration  in  a  way  which  will  be  described  in  the  following  section. 

The  calculations  that  have  been  described  in  this  section  are  effected  in 
the  subroutine  WF  BNDV  (or  in  the  entry  DECOMP  in  WF  BNDY). 


IV.  EXCITATION  FACTORS 


In  this  section  the  excitation  factors,  A,  which  have  been  programmed  are 
summarized  and  in  section  VII  their  usage  in  mode  sum  calculations  is 
discussed.  The  formulas  have  been  derived  in  reference  [8]  and  include 
"height  gain"  effects  associated  with  transmitter  and  receiver.  In  that  sense 
they  depart  from  the  conventional  definition.  The  formulas  programmed  differ 
from  those  of  reference  [8]  only  in  choice  of  sign  and  by  virtue  of  using  here 
the  complement  of  the  eigenangle  used  in  reference  [8].  The  formulas  are: 


1/2 

A1  =  fr790-e’1;l/4(1  +  "S")Hy(ZR)^1^R^^)hy+iR"^iey> 


Ao  =  - e1  37t/4( l+iR1  )E  (Zp)(,iRiiiRiih"VlR|lllR,1e1'") 

4  S4/  dF/30  y  R  y  y 


A3  =  - e"in/4(l+11R1()E  (2„)((l-1Riji)hVR11j1e;) 

s^aF/ae  x  K  y  y 


(22) 


In  the  above,  F  is  the  modal  function  (i.e.,  the  determinant  of  1  -  RR).  The 
superscript  m  denotes  the  source  while  the  subscripts  on  A  denote  the 
component  of  the  electric  field  at  the  receiver  for  which  the  excitation 
factor  applies.  The  convention  governing  m  and  j  is: 


Am 

J 


j  =  1  z  component  m  = 
j  =  2  y  component  m  = 
j  -  3  x  component  m  = 

m  = 
m  = 
m  = 


1  vertical  electric  dipole 

2  horizontal  electric  dipole  broadside 

3  horizontal  electric  dipole  end  fire 

4  vertical  magnetic  dipole 

5  horizontal  magnetic  dipole  broadside 

6  horizontal  magnetic  dipole  end  fire 


The  R  and  R  in  equation  (22)  have  been  defined  in  the  previous  section.  The 
quantities  h™  and  e™  are  the  full  wave  y  components  of  the  magnetic  and 
electric  fields  of  the  rf  wave  and  are  dependent  upon  the  type  and  location  of 


the  source  as  indicated  by  the  superscript  m.  The  dependence  on  altitude  of 

the  receiver,  z^,  in  the  input  coordinate  system  is  given  by  the  functions 

H  (zD),  Ew(z„)  and  E  (zD).  A  complete  description  of  these  functions  is  given 
y  k  y  k  x  k 

in  reference  [3].  Finally,  S  and  C  are  the  sine  and  cosine,  respectively,  of 
an  eigenangle,  6  ,  and  8F/30  is  the  derivative  of  the  mode  function  evaluated 
at  9  . 
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V.  PROGRAM  DESCRIPTION 


This  section  describes  the  subroutines  in  the  SATELLITE  program  listed  in 
appendix  A.  Many  of  the  subroutines  are  only  slight  modifications  of  those 
given  in  reference  [3].  The  subroutines  WKB,  WKBVAR,  QGAMMA  and  DDKXMT  have 
been  added  for  the  purpose  of  implementing  the  WKB  formalism.  Principal  out¬ 
put  of  the  program  is  the  excitation  factors  defined  in  section  IV.  A  chart 
showing  the  essential  structure  of  the  SATELLITE  program  follows  on  pages  22 
through  25. 

SUBROUTINE  MAIN 

MAIN  first  calls  for  input  of  ionic  species  data  in  XINPUT  and  for 
computation  of  height  gain  components  and  excitation  factors  via  WAVFLD. 
After  executing  WAVFLD,  excitation  factors  are  available  in  the  array  EFIELD 
(M,J)  where  M  is  the  source  designation  and  J  the  field  component  designation 
consistent  with  the  convention  given  in  section  IV. 

SUBROUTINE  XINPUT  (ISTART,  ISTOP) 

XINPUT  controls  read-in  of  input  parameters  via  NAMELIST  statements  and 
ionic  species  densities  and  collision  frequencies  as  a  function  of  altitude. 
Common  areas  are  set  up  as  required.  ISTART  is  set  to  1  before  the  first  call 
to  XINPUT  and  to  0  upon  subsequent  calls.  ISTART  =  1  implies  all  necessary 
data  are  to  be  read  in  and  ISTART  =  C  signals  that  previously  read  data  are  to 
be  updated,  with  all  unspecified  parameters  remaining  unchanged.  If  a  value 
ISTOP  =  0  is  returned  by  XINPUT,  then  more  input  data  are  specified  in  the 
data  deck  for  subsequent  calls  to  WAVFLD,  whereas  if  a  value  ISTOP  =  1  is 
returned,  the  data  read  were  the  last  data  in  the  data  deck,  so  that  XINPUT 
should  not  be  called  again.  XINPUT  effects  the  transformation  from  the  input 
coordinate  system  (x,  y,  z)  to  the  prime  coordinate  system  (x1  ,  y 1  ,  z1). 

The  data  deck  is  divided  into  several  parts,  each  of  which  is  marked  by 
an  identifier  card  with  the  identifier  DATUMFOL,  PROFILE,  COLLFREQ,  QUIT  or 
STOP.  Each  of  these  identifiers  is  described  in  the  following  section. 
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SUBROUTINE  WAVFLD  (EX,  EY,  EZ,  HX,  HY,  HZ) 

WAVFLD  calls  for  Runge-Kutta  integration  of  two  independent 
(corresponding  to  outgoing  wave  in  the  prime  system  at  TOPHT')  field  solutions 
from  TOPHT1  to  LWSTHT'  at  DELHT  increments.  The  two  solutions  are  combined  to 
satisfy  the  modal  polarization  condition  and  decomposed  into  magneto-ionic 
components  at  LWSTHT1.  At  LWSIHT'  the  full  wave  whistler  component  is  matched 
to  the  corresponding  WKB  component  and  the  solution  carried  to  WKBHT1  via  WKB 
formulas.  At  WKBHT'  the  strength  of  the  solution  is  determined  for  each 
source  and  orientation.  WAVFLD  then  calls  for  the  back  substitution  of 
normalizing  values  which  have  been  saved  in  WF  STOR.  The  excitation  factors 
are  then  calculated  in  WAVFLD  and  placed  in  the  array  EFIELD(M,J).  The  field 
components  placed  at  DELHT  intervals  in  the  arrays  EX,  EY,  EZ,  HX,  HY,  HZ  are 
extraneous  output  in  the  present  application. 

SUBROUTINE  WFINTG  (TOPHT1,  LWSTHT',  DELHT,  I  FLAG) 

WFINTG  performs  Runge-Kutta  integration  of  the  two  solution  vectors  in  P 
from  TOPHT'  to  LWSTHT1.  Numerical  solutions  are  saved  at  DELHT  increments. 
Accuracy  is  maintained  by  continually  adjusting  the  step  size  used  in  the 
numerical  integration.  The  current  step  size  (call  it  h)  is  used  to  obtain  an 
estimate  of  P,  and  then  a  better  estimate  is  obtained  by  using  two  integra¬ 
tions  with  step  size  h/2.  If  the  two  solutions  agree  within  an  error  of 
PRECSN  (an  input  parameter  normally  set  to  3.D-5),  the  better  estimate  is 
accepted.  The  step  size  is  automatically  decreased  to  h/2  if  the  two 
estimates  differ  by  more  than  PRECSN,  and  the  integration  is  repeated.  If  the 
error  is  significantly  greater  than  PRECSN,  a  step  size  h/2  is  used.  If  the 
error  is  significantly  less  than  PRECSN,  the  step  size  2h  is  used  if  it  also 
yields  an  error  less  than  PRECSN.  These  tests  thus  form  an  automatic  step 
size  correction.  In  this  program  the  call  to  WFINTG  sets  IFLAG  equal  to  zero. 

ENTRY  INIT  T 


INIT  T  is  an  entry  in  subroutine  TMTRX.  INIT  T  sets  up  height 
independent  values  to  be  used  in  T  matrix  calculation.  These  include  the 


internally  set  flag  ISO  (ISO  =  1  for  isotropic  calculation,  0  otherwise). 
Also  set  are  the  angular  radio  frequency,  the  wave  number,  direction  cosines 
of  the  geomagnetic  field,  the  complex  sine  and  cosine  of  THETA. 


SUBROUTINE  TMTRX(HT') 

TMTRX  computes  the  value  of  the  T  matrix  and  a  specific  height  HT '  (as 
usual  the  prime  simply  implies  that  HT  is  referred  to  the  prime  coordinate 
system).  The  T  matrix  depends  upon  input  ionospheric  parameters  (species 
density,  collision  frequency,  angle  of  propagation,  magnetic  field,  etc.). 
The  susceptibility  matrix,  M,  for  each  species  in  the  ionosphere  is  computed, 
the  effect  of  earth  curvature  is  included  and  the  T  matrix  is  computed  from 
the  susceptibility  matrix  elements.  The  equations  used  to  evaluate  M  and  T 
are  given  in  Clemmow  and  Heading  [6], 

SUBROUTINE  WFDENS  (HT1 ,  EN,  COLL) 

WFOENS  computes  the  species  density  (returned  in  EN)  and  collision 
frequency  (returned  in  COLL)  at  height  HT'  for  up  to  five  species  in  the 
ionosphere.  EN  and  COLL  are  determined  from  corresponding  profiles  in  the 
common  field  WFPROF.  LHT  and  MHT  are  integer  values  which  indicate  which 
profile  values  are  to  be  interpolated  to  find  values  at  the  height  HT*.  The 
EN  (or  COLL)  values  are  given  by  logarithmic  interpolation  of  the  profile 
values  of  heights  ENHF(LHT  +  1)  and  ENHT(LHT)  or  (COLLHT(MHT  +  1)  and 

COLLHT(MHT)). 

SUBROUTINE  WF  INIT(P) 

WF  INI!  sets  the  two  starting  field  vectors  P,  at  TOPHT'  subject  to  the 
condition  of  outgoing  wave  in  the  prime  system.  The  two  starting  solutions 
correspond  to  TE  and  TM  polarizations. 

SUBROUTINE  PUER IV( P , DPDH ) 

PDERIV  computes  the  height  derivatives  of  the  two  field  vectors  in  P 
according  to  Clemmow  and  Heading  [6].  The  derivative  is  returned  in  DPDH. 
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SUBROUTINE  XFER  (A,  B,  M) 


Transfers  the  N  element  array  A  into  B. 

ENTRY  WF  STOR 

This  is  an  entry  in  WFSCAL  where  certain  values  obtained  during 
integration  through  the  ionosphere  are  saved  for  later  use.  The  solution 
matrix  P,  the  height  for  which  P  is  a  solution  and  a  height  integer  index  are 
saved  along  with  orthogonal i zation  and  normalization  values  OSUM,  APROD  and 
BPROD.  In  addition,  values  of  the  susceptibility  matrix  elements  M31,  M32  and 
M33,  which  are  needed  to  compute  the  EZ  component  of  the  electric  field,  are 
saved  at  each  height. 

SUBROUTINE  WF  STEP  (P.DPDH,  HT ' ,  DELHT ,  IFLAG) 

WF  STEP  numerically  advances  the  solution  matrix  P,  using  the  derivative 
DPDH,  from  HT'  to  HT1 -DELHT  by  the  Runge-Kutta  method.  IFLAG,  set  internally, 
controls  the  calculations  at  intermediate  points  between  HT'  and  HT'-DELHT  at 
which  evaluations  are  required  for  comparison  of  the  second  and  fourth  order 
Runge-Kutta  integrations. 

SUBROUTINE  WF  SCAL  (PP,  IFLAG) 

WF  SCAL  normalizes  and  orthogonal izes  the  solution  vectors  PP  according 
to  the  formulas  of  reference  [4],  This  scaling  must  later  be  removed  to  yield 
correct  unsealed  solutions.  Control  for  calculating  the  quantities  needed  for 
removal  of  the  scaling  is  the  internally  set  IFLAG.  Calling  WF  SCAL  sets 
IFLAG  =  0  in  this  program. 

ENTRY  DECOMP(HT') 

This  is  an  entry  in  WF  BNDY.  The  two  full  wave  solution  vectors,  P,  are 
decomposed  at  HT1  =  LWSTHT'  into  upgoing  (in  the  prime  system)  magneto-ionic 
components  and  combined  to  satisfy  the  modal  polarization  condition  given  in 
equation  (19). 
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SUBROUTINE  EIGVAL  (A,  SIG) 

EIGVAL  computes  the  eigenvalues  (returned  in  SIG)  or  the  arbitrary  4*4 
complex  matrix  A.  The  characteristic  polynomial  P(A)  of  A  is  determined  by 
explicitly  computing  det(A-AI)  =  P(A).  The  roots  of  this  quartic  polynomial 
(as  computed  in  closed  form,  see  subroutine  QUARTO)  are  the  desired  eigen¬ 
values  of  A.  If  the  ionosphere  is  isotropic,  the  eigenvalues  of  A  are 
computed  as  the  roots  of  two  quadratics. 

SUBROUTINE  QUARTO  (FOUR  B3,  SIX  B2,  FOUR  Bl,  B0,Q) 

QUARTO  computes  the  four  roots  of  the  polynomial  Q4  +  FOUR  B3  *  Q3  + 
SIX  B2  *  Q2  +  FOUR  Bl  *  Q  +  BO  in  closed  form  [9].  Up  to  ten  applications  of 
Newton's  iteration  are  then  performed  to  improve  ihe  accuracy  of  the  roots,  if 
necessary. 

SUBROUTINE  EIGVEC  (A,  SIG,  VEC) 

EIGVEC  computes  the  four  eigenvectors  of  A  (returned  in  the  four  columns 
of  the  matrix  VEC).  The  eigenvalues  of  A  are  assumed  given  in  SIG.  Each 
eigenvector  X.  is  computed  as  the  solution  of  (A  -  o.I^.  =  0.  Since  det(A  - 
(i.  I)  =  0,  this  system  has  a  solution  if  one  element  of  X^  is  chosen 
arbitrarily.  In  order  to  reduce  the  special  isotropic  case,  the  first  element 
of  X ^  and  X^  are  arbitrarily  set  to  A^-  The  other  three  elements  of  each 
vector  are  then  determined  as  the  solution  of  a  linear  system  in  three 
unknowns.  If  the  ionosphere  is  isotropic,  the  eigenvectors  (corresponding  to 
the  quadratic  roots  described  previously)  are  computed  in  simplified  form. 

SUBROUTINE  WF  SORT  (Q,  A  IFAIL) 

WF  bORT  arranges  the  eigenvalues  Q  and  eigenvectors  A  of  the  Booker 
quartic  so  that  they  occur  in  the  order  of  upgoing  fast  (evanescent),  upgoing 
slow  (travelling  or  whistler  component),  downgoing  slow,  and  uowngoing  fast  in 
the  prime  coordinate  system.  If  the  eigenvalues  are  too  large  in  magnitude  or 
if  they  cannot  be  ordered  IFAIL  is  set  to  1  and  an  error  message  results. 


SUBROUTINE  CIMVER  (A,  AINV,  N,  NDIM,  ERR) 

CINVER  computes  the  inverse  (returned  in  AINV)  of  the  NxN  matrix  A.  NDIM 
is  an  integer  variable  which  must  be  greater  than  or  equal  to  N.  ERR  is  the 
estimated  relative  error  of  the  inverse  matrix. 

SUBROUTINE  CLINEQ  (A,  B,  X,  N,  NDIM,  I  FLAG,  ERR) 

CLINEQ  computes  the  solution  of  simultaneous  linear  equations  with 
complex  coefficients.  That  is,  it  solves  the  matrix  A  *  X  =  B  for  the  vector 
X  of  length  N,  given  the  matrix  A  of  size  N  by  N  and  the  vector  B  of  length  N. 
The  matrix  A  is  destroyed  by  CLINEQ.  NDIM  is  an  integer  variable  which  must 
be  greater  than  or  equal  to  N.  IFLAG  is  an  integer  variable  normally  set  to 
0.  ERR  indicates  the  relative  error  in  the  computed  solution  of  vector  X. 

SUBROUTINE  RBARS  (C,  S,  RBAR11 ,  RBAR22) 

RBARS  calculates  the  plane  wave  reflection  coefficients  „RM=RBAR11  and 
1Ri.=RBAR22  looking  towards  the  ground  from  LWSTHT  in  the  input  system. 
Evaluations  are  in  terms  of  modified  Hankel  functions  of  order  1/3.  The 
cosine  (C)  and  sine  (S)  of  the  eigenvalue,  THETA,  are  used  in  the  R 
determination. 

SUBROUTINE  MDHNKL  (Z,  HI,  H2,  H1PRME ,  H2PRME) 

MDHNKL  calculates  for  argument  Z  two  independent  solutions  (HI  and  H2) 
and  their  derivatives  (H1PRME,  H2PRME)  of  Stokes'  equation  by  methods 
described  in  reference  [10].  HI  and  H2  are  the  modified  Hankel  functions  of 
order  1/3  referred  to  above. 

SUBROUTINE  WKB 

If  STPFLG.NE.l,  WKB  extends  calcuation  of  the  solution  vector  from 
LWSTHT'  to  WKBHT'  by  means  of  the  formulas  in  section  II.  Integration  of  the 
phase  factors  q^  and  F^  (j  index  for  upgoing,  in  prime  system,  whistler  mode) 
is  performed  by  Simpson  rule  routine  with  fixed  step  size  that  can  be  either 
input  to  the  program  or  set  internally. 
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SUBROUTINE  WKVAR 

If  STPFLG.EQ.l,  WKBVAR  extends  calculation  of  the  solution  vector  from 

LWSHT'  to  WKBHT'  via  the  WKB  formulas  of  section  II.  Although  generally 

requiring  more  CAU  time  this  is  the  preferred  option  since  it  involves  checks 

on  step  size.  Integration  of  the  phase  factors  q.  and  r..  (j  index  for 

J  J  J 

upgoing,  in  prime  system,  whistler  mode)  is  performed  by  a  Simpson  rule 
routine  which  maintains  precision  by  adjustment  of  the  step  size  much  like  the 
Runge-Kutta  step  size  adjustment  in  WFINTG. 

ENTRY  INITDT 

This  is  an  entry  in  DDKXMT  where  all  height  independent  quantities  are 
initialized  before  extending  field  calculations  from  LWSTHT'  to  WKBHT'  via  the 
WKB  method. 

SUBROUTINE  Q  GAMMA  (HT\  DELHT ,  LWSTHT1,  WKBHT1,  Q,  GAMMA) 

Q  GAMMA  determines  the  Booker  quartic  solutions  for  the  least  attenuated 
upgoing  (in  the  prime  system)  magneto- ioni c  component  at  LWSTHT1  and  computes 
the  a.'s,  b.'s,  al's,  bl's,  A.,  F.  and  T- .  according  to  formulas  of  section 
II.  Full-wave  solutions  are  matched  at  LWSTHT'  to  the  WKB  solutions.  At 
other  heights  and  up  to  WKBHT'  coefficients  of  the  exponentials  in  equation 
(9)  of  section  II  are  calculated  along  with  coefficients  required  for 
calculating  the  EZ  and  HZ  fields. 

SUBROUTINE  DDKXMT 

DDKXMT  calculatzs  susceptibility  matrix  elements,  M,  the  T  matrix 
elements  and  their  derivatives  with  respect  to  height  in  the  height  range 


LWSTHT ' 


WKBHT' . 


WF  BNDY(B) 

This  computes  coefficients  BP  (I,2)--see  equation  (21)--for  each  source 
and  orientation  at  WKBHT'.  The  coefficients  are  used  to  determine  the 
strength  of  the  full  wave  solution  which  has  been  carried  from  within  the 
guide  to  WKBHT'  by  the  combination  of  full  wave  and  WKB  methods. 


ENTRY  WF  LOAD 


WF  LOAD  is  an  entry  in  WF  SCAL  which  makes  available  for  the  purpose  of 
back  substitution  into  the  full  wave  solutions  the  orthogonal izi ng  and 
normalizing  coefficients  saved  in  WF  STOR. 

ENTRY  HT  GAIN 

This  is  an  entry  in  RBARS  where  the  quantities  required  to  calculate  the 
properly  normalized  height  gains  at  the  receiver  altitude  are  evaluated. 


MAIN 


1 


Calls  for  input  data  via  XINPUT. 
Calls  for  computation  by  WAVFLD  of 
excitation  factors  for  the  ELF/VLF 
electric  field  components  beneath 
the  ionosphere  produced  by  elec¬ 
tric  and  magnetic  dipoles  at 
sate  1 1 i te  al ti tudes. 

XINPUT 

Input  electron  and  ion 
species  densities  and  col¬ 
lision  frequencies  as  a 
function  of  altitude  in  the 
input  coordinate  system  (z 
positive  into  the  iono¬ 
sphere).  Converts  input  to 
the  primed  coordinate 
system  (z1  positive  towards 
the  ground). 


WAVFLD 

Calls  for  Runge-Kutta  integration 
of  two  independent  field  solutions  from 
TOPHT'  to  LWSTHT'  at  DELHT  increments. 
The  two  solutions  are  combined  to 
satisfy  the  modal  polarization  condi¬ 
tion  and  decomposed  into  magneto- i oni c 
components  at  LWSTHT'.  WAVFLD  calls 
for  the  back  substitution  of  ortho¬ 
gonal  i  zing  and  normalizing  values 
(saved  as  data  in  WFSTOR).  It  calls 
for  the  extension  of  the  fields  from 
LWSTHT'  to  WKBHT '  via  Budden's  WKB 
formalism.  At  WKBHT'  the  field  com¬ 
ponents  are  matched  to  the  cor¬ 
responding  components  generated  by  each 
source.  Excitation  factors  for  the 
ELF/VLF  electric  field  components  are 
then  computed  for  an  input  receiver 
altitude  which  must  be  within  the 
isotropic  space  between  the  ground  and 
ionosphere. 
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WKB  or  WKBVAR 


INIT  DT 

This  is  an  entry 
in  DDKXMT.  All 
height  independent 
quantities  are  ini¬ 
tialized  before 
extending  field  cal¬ 
culations  from 
LWSTHT*  to  WKBHT’  via 
Budden's  WKB  formulas 
for  an  anisotropic 
medium. 


Q  GAMMA 

Beginning  at  LWSTHT1  in  the  prime 
system  computes  a  -'s,  b  .'s,  A.,  F  -,  r  - . , 

J  J  J  J  J  J 

aVs,  b'-'s  Cj  is  index  for  least 

J  J 

attenuated  outgoing  wave  in  the  prime 
system)  according  to  the  formulas  of 
section  II.  Matches  full  wave  fields 
of  LWSTHT'  to  WKB  solutions. 


DDKXMT 

Calculates  susceptibility  matrix 
elements,  M,  the  T  matrix  elements 
and  the  derivative  of  the  M,  T 
elements  with  respect  to  height. 


QUARTC 


WF  DENS 


VI.  SAMPLE  INPUT  AND  OUTPUT 


A.  INPUT 

Altitude  independent  parameters,  species  densities  and  collision 
frequencies  as  a  function  of  altitude  are  supplied  via  an  input  data  deck  on  a 
standard  input  unit.  Read-in  occurs  in  the  subroutine  XINPUT.  The  data  deck 
is  divided  into  several  parts,  each  of  which  is  marked  by  control  cards 
DATUMFOL,  PROFILE,  COLLFREQ,  QUIT,  or  STOP.  Each  of  these  control  cards  is 
described  below: 


(1)  DATUMFOL 

DATUMFOL  -  is  a  control  card  signifying  that  input  for  the 
namelist  DATUM  follows. 

&DATUM  -  initiates  reading  of  input  data  which  do  not  vary 
with  altitude.  The  data  is  read  in  namelist  input 
format.  All  altitude  related  quantities  are 
referenced  to  the  input  coordinate  system  (z 
positive  into  the  ionosphere,  ground  at  z  =  0). 
These  data  are: 


THETA  -  complex  angle  of  incidence  in  degrees  as  measured  at 
height  H.  THETA  must  be  applied  from  a  waveguide 
program  (e.g.,  ref.  [5]). 

FREQ  -  frequency  in  kHz. 

TOPHT  -  height  to  which  full  wave  integrations  are  carried 
(km). 

-  beginning  height  for  full  wave  integrations  (km). 
Also  equivalent  to  level  D  where  mode  equation  is 
eva 1 uated. 


?b 


LWSIHT 


tfKBHT 


lower  boundary  of  homogeneous,  anisotropic  slab  in 
which  the  transmitter  is  imbedded  (km).  Also  height 
to  which  field  solution  is  propagated  via  WKB 
formalism  from  TOPHT. 


DELHT  -  height  increment  in  km's  at  which  field  solutions 
are  saved. 

PRECSN  -  accuracy  to  be  maintained  locally  in  the  numerical 
integrations.  Usually  taken  to  be  the  default  value 
of  3.0E-5. 

AZIM  -  azimuth  of  propagation  path  in  degrees,  measured 
clockwise  from  geomagnetic  north. 

CODIP  -  codip  of  geomagnetic  field  in  degrees. 

MAGFLD  -  geomagnetic  field  strength  in  webers/square  meter. 

C0EFNU(5)-  coefficient  of  exponential  form  of  collision 
frequency  (if  not  specified  by  profile).  Up  to  five 
values  may  be  specified.  One  for  each  species. 

EXPNU(5)  -  exponent  of  exponential  form  of  collision  frequency 
(if  not  specified  by  profile).  Up  to  five  values 
may  be  specified.  One  for  each  species. 

ALPHA  -  earth  curvature  coefficient  in  inverse  km.  Default 
is  3.14E-4. 

SIGMA  -  ground  conductivity  in  siemens/meter.  Default  value 
is  seawater  value  of  4.64. 


EPSLON  -  ground  permittivity  in  farads/meter.  Default  value 
is  7. 172015D-10.  This  corresponds  to  a  dielectric 
constant  of  81  for  seawater. 


EPSLNO  -  permittivity  of  free  space  (8.8544E-12)  in 
farads/meter. 

H  -  altitude  in  km  at  which  modified  refractive  index  is 

unity.  Eigenangles  are  referenced  to  this  altitude. 
H  must  be  consistent  with  the  H  setting  in  the  wave¬ 
guide  program  which  supplied  the  eigenangles. 

TXHT  -  transmitter  altitude  (km). 

RXHT  -  receiver  altitude  (km). 

ITR  -  integer  flag  which  should  be  set  to  zero  in  this 

program. 

RMAG(4)  -  magnitude  of  the  plane  wave  reflection  elements  of 
R.  R ( 1 )  =  1(R„,  R(2)  =  ^Rj. ,  R(3)  =  j.R„  ,  R(4)  =  „RA. 

These  elements  are  referenced  to  level  D  (or 

equivalently  LWSTHT)  in  the  input  coordinate  system 
and  must  be  supplied  from  a  waveguide  program  such 
as  that  of  reference  [5]. 

RANG(4)  -  phase  of  the  plane  wave  reflection  elements  of  R. 

These  elements  are  referenced  to  level  D  (or 

equivalently  LWSTHT)  in  the  input  coordinate  system 
and  must  be  supplied  from  a  waveguide  program  such 
as  that  of  reference  [5]. 

XTRMAG  -  magnitude  of  the  ground  based  vertical  dipole 

excitation  factor  for  the  vertical  electric  field, 
Ez.  This,  along  with  XTRANG,  is  used  for  the 

purpose  of  obtaining  the  derivative  of  the  modal 
equation,  dF/dO.  XTRMAG  must  be  supplied  from  a 
waveguide  program  such  as  that  of  reference  [Sj. 


XTRANG  -  phase  of  the  ground  based  vertical  dipole  excitation 
factor  for  the  vertical  electric  field,  Ez.  This, 
along  with  XTRMAG,  is  used  for  the  purpose  of 
obtaining  the  derivative  of  the  modal  equation, 
dF/de.  XTRANG  must  be  supplied  from  a  waveguide 
program  such  as  that  of  reference  [5]. 

NUMDIV  -  number  of  divisions  into  which  DELHT  is  divided  when 
using  subroutine  WKB  (i.e.  ,  when  STPFLG  ^  1)  with 
fixed  input  step  size.  NUMDIV  must  then  be  a 
multiple  of  2.  If  STPFLG  i-  1  and  NUMDIV  =  0  then 
the  program  will  determine  a  constant  step  size  to 
be  used  with  the  Simpson  rule  integration.  It  is 
repeated  here  that  the  preferred  option  is  with 
STPFLG  =  1  because  of  the  checks  on  the  adequacy  of 
the  step  size. 

STPFLG  -  integer  flag  which  must  be  set  to  1  if  the 
recommended  WKBVAR  subroutine  is  to  be  used  to 
effect  the  WKB  integrations.  Other  options  are 

discussed  above  (NUMDIV). 

NRSPEC  -  number  (integer)  of  species  in  the  ionosphere.  Can 
take  on  values  up  to  5.  Default  value  is  1. 

CHARGE(5)-  sign  of  charge  of  each  species  in  proton  units.  For 
an  electron,  the  charge  is  -1.0.  Default  values  are 
(-1.0,  1.0,  -1.0,  1.0,  -1.0). 

RATI0M(5)-  mass  of  each  species  relative  to  mass  of  an 
electron.  Default  values  are  (1.0,  5.8D4,  5.8D4, 
5.8D4,  5.8D4). 

&END  -  signifies  end  of  the  DATUM  namelist  input. 


PROFILE 


PROFILE  -  is  a  control  card  initiating  reading  of  the 
ionospheric  profile  cards.  All  altitudes  and 
related  quantities  are  referenced  to  the  input 
coordinate  system  (z  positive  into  the  ionosphere, 
ground  at  z  =  0).  The  control  card  PROFILE  is 
followed  by  an  alpha-numeric  card  which  is  used  to 
identify  the  profile.  The  profile  is  input  starting 
at  the  top  of  the  ionosphere  (i.e.,  WKBHT).  The 
cards  must  be  input  in  descending  order  by  height. 
The  profile  cards  contain  the  height  in  kilometers 
and  the  species  densities  in  particles  per  cubic 
centimeter.  A  maximum  of  five  species  can  be 

specified.  In  the  special  case  of  three  species, 
only  two  are  specified  on  the  card.  The  first  is 
assumed  to  be  electrons  and  the  second  is  positive 
ions.  The  third  species,  negative  ions,  is 
calculated  by  subtracting  the  electron  density  from 
the  positive  ion  density.  All  three  species  are 
listed  on  the  computer  printout.  If  the  value  of 
any  species  density  is  less  than  or  equal  to  zero, 
it  is  set  in  the  program  to  1.0E-40.  The  heights 
are  punched  in  the  form  xxx.xx  with  the  decimal 
point  in  column  5,  and  the  densities  are  punched  in 
the  form  x.xxD+xx  with  the  decimal  points  in  column 
15,  25,  etc.  All  species  data  must  be  specified 
except  for  the  special  case  discussed  above.  The 
end  of  the  profile  is  indicated  by  a  dummy  height  of 
359.99. 
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(iii)  COLLFREQ 


COLLFREQ  -  is  a  control  card  initiating  reading  of  the 
collision  frequency  profile.  All  altitudes  and 
related  quantities  are  referenced  to  the  input 
coordinate  system  (z  positive  into  the  ionosphere, 
ground  at  z  =  0).  The  control  card  COLLFREQ  allows 
for  using  nonexponential  collision  frequencies.  A 
strictly  exponential  collision  frequency  may  be 
specified  in  namelist  input  by  the  variables  COEFNU 
and  EXPNU.  If  a  collision  frequency  profile  deck  is 
used  it  overrides  COEFNU  and  EXPNU.  Collision 
frequencies  for  all  species  must  be  input.  The  card 
preparation  is  just  as  above  with  the  species 
density  values  replaced  by  collision  frequency 
values  (in  collisions/sec)  on  the  cards.  The  end  of 
the  profile  is  indicated  by  a  dummy  height  of 
999.99. 


(iv)  QUIT 

QUIT  -  is  a  control  card  which  indicates  the  end  of  input 
data  for  a  call  to  WAVFLD.  After  calling  WAVFLD, 
XINPUT  can  be  read  and  used  as  data  for  a  subsequent 
call  to  WAVFLD.  This  allows  several  runs  to  be  made 
with  one  input  deck.  Note  that  only  those  data 
which  are  changed  need  be  specified  after  the  card 
QUIT. 

(v)  STOP 

STOP  -  control  card  which  indicates  that  there  are  no  more 
input  data  and  that  the  run  is  to  be  terminated 
after  the  next  call  to  WAVFLD.  Note  that  if  QUIT  is 
encountered  ISTOP  is  set  to  zero,  but  if  STOP  is 
encountered  ISTOP  is  set  to  one. 


A  schematic  of  an  input  deck  is  shown  on  page  34,  and  one  actual  sample  input 
is  shown  on  pages  35  through  37. 

B.  OUTPUT 

The  output  shown  on  pages  38  through  41  corresponds  to  the  first  DATUMFOL 
shown  on  page  36.  The  output  corresponds  to  the  STPFLG  =  1  setting.  Output 
corresponding  to  the  additional  DATUMFOL  have  not  been  shown  simply  to  cut 
down  on  the  number  of  pages. 

First  come  NAMELIST  and  profile  printouts  in  the  input  coordinate  system 
(z  positive  into  the  ionosphere,  ground  at  z  -  0).  Natural  logarithms  of  the 
profiles  in  the  prime  system  are  also  printed.  Next  follows  output  which 
gives  the  number  of  Runge-Kutta  integration  steps  used  in  WAVFLD  during  the 
course  of  integrating  from  TOPHT '  to  LWSTHT '  as  well  as  the  smallest  and  the 
average  step  sizes  used  in  the  WKB  extension  from  LWSTHT'  to  WKBHT'.  This  is 
followed  by  the  principal  output;  two  nine  element  arrays  which  are  the 
excitation  factors  for  both  electric  and  magnetic  dipole  transmitters.  The 
column  headings  indicate  excitation  factors  for  the  z,  y  and  x  electric  field 
components  in  the  input  coordinate  system.  The  row  labels  indicate  excitation 
by  vertical  (V),  horizontal  broadside  (HB)  and  horizontal  end-on  (HE)  point 
dipoles  of  either  electric  or  magnetic  type. 

We  list  below  a  comparison  at  20  kHz  between  the  excitation  factors 
obtained  using  full  wave  (FW)  methods  up  to  500  km  as  opposed  to  using  WKB 
formalism  between  125  km  and  500  km.  The  agreement  will  be  seen  to  be 
exce 1  I ent. 

ELECTRIC  DIPOLE  TRANSMITTER 


ANG  (RAD.  ) 


ANG  (RAD. ) 


MAG  ANG  (RAD.) 


EXAMPLE  OF  INPUT  DECK 

DATUMF0L 

&DATUM 
( i nput  data) 

&END 


PR0FILE 

(profile  cards  for  each  specie) 

999.99 

C0LLFREQ 

(coll,  frequency  cat'ds  for  each  specie) 

999.99 

QUIT  (end  of  input  for  this  run) 

DATUMF0L 

&DATUM 

(changes  in  input  data) 

&END 

PR0FILE 

(specify  entire  new  profile  deck) 

999.99 

QUIT  (end  of  input  for  this  run) 

DA1UMF0L 

&DATUM 

(changes  in  input  data) 

At  NO 

3I0P  (end  of  data  deck) 
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V.  MODE  SUMS 

The  excitation  factors  defined  in  the  previous  section  are  useful  for 
calculating  the  strengths  of  the  electric  field  components  in  the  earth- 
ionosphere  guide.  In  this  section  explicit  expressions  are  given  for 
calculating  those  field  components.  The  formulas  given  are  asymptotic  and 
restricted  to  the  range 

|kpsin@^|>>l  and  | k(na-p)si n0^| >>1 

with  k  the  free  space  wave  number,  0^  the  eigenangle  for  Mode  N,  a  the  earth's 
radius  and  p  the  transmitter  receiver  great  circle  separation.  The  mode  sums 


C,I£f: 


rm/pv\_  1  kHz  r  ,m  _-ikpsin0„ 

v  i :  i:i:l 


Emftjv 
j  \m 


C  IAf~^ 
L21ArkHz 


^  ■>4Z,._  1  a”  e'ikpSin6N 

sin/Eji172  N  jN 
\  a  /  / 


j  =  1,2,3 
m  =  4,5,6 


where  the  sum  is  over  the  number  of  modes  and 


Cj  -  2.86  x  10 


C2  =  5.98  x  10 


f^  =  frequency  in  kHz 


l£  =  electric  dipole  current  moment  -  Amp-m 


IA  =  magnetic  dipule  current  moment  -  Amp-m2 


A.  -  excitation  factors  defined  in  previous  section. 
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The  field  strengths  E.  are  in  terms  of  microvolts  per  meter.  A  program  for 

J 

computing  and  plotting  the  mode  sums  is  included  in  appendix  B.  The  latter 
takes  punched  output  from  the  excitation  factor  program  given  in  appendix  A. 
Results  of  the  mode  sum  program  are  shown  in  figures  1  through  3  for  electric 
dipole  (ED)  and  magnetic  dipole  (MD)  excitation  of  EX,  EY  and  EZ.  The  dipole 
moments  in  this  example  have  been  chosen  to  correspond  to  one  kilowatt  free 


;iRWR 1 1  TO  SRN  D'EGO  GLOBRL  NIGHT  PROFILE 

PREQ-20.00Q  KHZ  TXHT-500 . 9  KM  RXHT-10.0  KM 
RZIM-53 . 5  DEGREES  CODIp-  39.0  DEGREES 
MRGplD-4.31  *10~5  rt/M‘  SIGMR-  4.64 ^lO'S'M  EpSR  =  31. 9 
IL”  i  .  69  *  10  *  R_M  IR- 4.Q3  <107R-M2 


DB/^n  CZ  -  DB/mV/M  BZ  -  DB/pV/n 
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SPRINGFIELD  VA  22150 
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US  ARMY  RESEARCH  O1-  p  ICE 
PO  BOX  122  1 1 

TRIANGLE  PARK,  NC  27/09 
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DEPARTMENT  OF  THE  NAVY 
CHIEF  Qp  NAVAL  OPERA  I  IONS 
NAVY  DEPARTMENT 
WASHINGTON.  DC  203  50 
OP  941 
OP-604C3 

OP943  (LCDR  HUFF) 

OP  981 

CHIEF  OF  NAVAL  RESEARCH 
NAVY  DEPARTMENT 
ARLINGTON.  VA  22217 
CODE  402 
CODE  420 
CODE  42) 

CODE  461 
CODE  464 

COMMANDING  OFFICER 
NAVAL  INTELLIGENCE  SUPPORT  CENTER 
4301  SUITLAND  RD  bLDG  5 
WASHINGTON.  DC  20390 

COMMANDER 

NAVAL  OCEAN  SYSTEMS  CENTER 
SAN  DIEGO  CA  92152 

CODE  532  !JH  RICHTER) 

CODE  532  IRICHARD  A  P  APPART- 
CODE  532  JERRY  F  FERGUSON) 

CODE  6224  IlR  hiTNEY: 

COMMANDING  OFFICER 
NAVA-  RESEARCH  LABORATORY 
WASHINGTON.  DC  20375 
CODE  54  10  (JOHN  DAVIS; 

CODE  7701  (JACK  D  SHOWN) 

CODE  5461  TRANS  IONOP«OP 
COOE  5465  PROP  APPLICATIONS 
CODE  5460  ELECTROMAG  PROP  BR 
CODE  2600  TECH  LIBRARY  (2 1 
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NAVY  UNDERWATER  SOUND  LABORATORY 
F OR T THUMPul  L 
NEW  LONDON.  CT  06320 
PETER  BANNISTER 
DA  Mil  LFR 

OIHE;  MOP 

STRATEGIC  SYSTEMS  PROJECT  OFFICE 
NAVY  DEPAR  rMENT 
W  AS:-.  iNGT ON  DC  20376 
MSP  -2141 

D i;P AH TMEN1  OF  the  AIR  =  ORCE 
COMMANDER 
A  DC  2)C 

FNT  AFB.  CO  80912 
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COMMANDER 
A DC  OM/XPD 
ENT  AFB.  CO  80912 
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XP 

AF  GEOPHYSICS  LABORATORY.  AFSC 
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CRU  (S  HOROWITZ) 

PHP  (JULES  AARONS) 

OPR  (JAMES  C  ULWICKI 
OPR  (ALVA  T  STAIR) 

SUOL  (RESEARCH  LIBRARY)  (2 

AF  WEAPONS  LABORATORY.  AFSC 
KlRTLAND  AFB,  NM  87117 

Sul  (? 

SAS  (JOHN  M  KAMVI) 

OYC  ICAPT  L  WITTWER) 

AFT  AC 

DATRICK  AFB.  FL  32925 
TN 
TO  3 
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WRIGHT  PATTERSON  APR.  Cih  .16433 

A  AD 

COMMANDER 

FOREIGN  TECHNOLOGY  DIVISION  AFSC 
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COMMANDER 

ROME  AIR  DEVELOPMENT  CENTER.  AFSC 
GRIEFISS  AFB,  NV  13440 
EMTLD  DOC  LIBRARY 

COMMANDER 

ROME  AIR  DEVELOPMENT  CENTER.  AFSC 
HANSCOM  AFB.  MA  01731 
EEP  JOHN  RASMUSSEN 

SAMSO/MN 

NORTON  AFB,  CA  92409 

MINUTEMAN  (NMMl  LTC  KENNEDY' 

COMMANDER  IN  CHIEF 
STRATEGIC  AIR  COMMAND 
OFFUTT  AFB.  NB  68113 
NRT 

XPFS  (MAJ  BRIAN  G  STEPHAN) 

DOK  (CHIEF  SCIENTIST) 

US  ENERGY  RESEARCH  AND  DEV  ADMIN 
DEPARTMENT  OF  ENERGY 
ALBUQUERQUE  OPERATIONS  OF FF ICE 
PO  BOX  6400 

ALBUQUERQUE.  NM  87115 

DOC  CON  FOR  DW  SHERWOOD 

division  of  military  application 

DEPARTMENT  OF  ENERGY 
WASHINGTON.  DC  20545 

DOC  CON  FOR  DONALD  I  GALE 

LAWRENCE  LIVERMORE  NATIONAL 

LABORATORY 

PO  BOX  808 

LIVERMORE.  CA  94550 
GLENN  C  WERTH  L-216 
TECH  INFO  DEPT  L-3 

LOS  ALAMOS  NATIONAL  SCIENTIFIC 

LABORATORY 

PO  BOX  1663 

LOS  ALAMOS.  NM  37545 

DOC  CON  FOR  TF  TASCHEK 
DOC  CON  FOR  DR  WESTERVELT 
DOC  CON  FOR  Pw  KEATON 
DOC  CON  FOR  JH  COON 

SANDIA  LABORATORY 

LIVERMORE  NATIONAL  LABORATORY 

PO  30X  969 

LIVERMORE.  CA  94550 

DOC  CON  FOR  BE  MURPHEY 
DOC  CON  FOR  TB  COOK  ORG  8000 
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SANDIA  NATIONAL  LABORATORY 
PO  BOX  6800 

ALBUQUERQUE.  NM  87115 

DOC  CON  FOR  SPACE  PROJ  DIV 

DOC  CON  FOR  AD  THORNBROUGH.  ORG  1245 

UOC  CON  FOR  WC  MYRA 

DOC  CON  FOR  314  1  SANDIA  RPT  COLL 

OTHER  GOVERNMENT 
DEPARTMENT  OF  COMMERCE 
NATIONAL  BUREAU  OF  STANDARDS 
WASHINGTON,  DC  20234 
RAYMOND  T  MOORE 

DEPARTMENT  OF  COMMERCE 
OFFICE  OF  TELECOMMUNICATIONS 
INSTITUTE  FOR  TELCOM  SCIENCE 
BOULDER,  CO  80302 
WILLIAM  F  UTLAUT 
LA  BERRY 
A  GLENN  JEAN 
DD  CROMBIE 
JR  WAIT 

DEPARTMENT  OF  TRANSPORTATION 
OFFICE  OF  THE  SECRETARY 
T AD-44.1 ,  ROOM  10402-B 
400  7TH  STREET,  SW 
WASHINGTON,  DC  20590 
RL  LEWIS 
RH  DOHERTY 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 
AEROSPACE  CORPORATION 
PO  BOX  92957 
LOS  ANGELES,  CA  90009 
IRVING  M  GARFUNKEL 

ANALYTICAL  SYSTEMS  ENGINEERING  CORP 
5  OLD  CONCORD  RD 
BURLINGTON,  MA  01803 
RADIO  SCIENCES 

THE  BOEING  COMPANY 
PO  BOX  3707 
SEATTLE,  WA  98124 
GLENN  A  HALL 
JF  KENNEY 

CUMPUTFR  SCIENCES  CORPORATION 
PO  BOX  530 

6565  ARLINGTON  BLVD 
FALLS  CHURCH.  VA  22046 
D  BLUMBERG 

UNIVERSITY  OF  DENVER 
COLORADO  SEMINARY 
DENVER  RESEARCH  INTSITUTE 
PO  BOX  10127 
DENVER.  CO  80210 
DONALD  DUBBERT 
HERBERT  REND 
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ESL .  INC 
495  JAVA  DRIVE 
Sunnyvale,  ca  94086 
JAMES  MARSHALL 

GENERAL  ELECTRIC  COMPANY 
SPACE  OIVISION 
VALLEY  FORGE  SPACE  CENTER 
GODDARD  BLVD  KING  OF  PRUSSIA 
PO  BOX  8555 

PHILADELPHIA,  PA  19101 

SPACE  SCIENCE  LAB  (MH  BORTNERi 

GENERAL  ELECTRIC  COMPANY 
TFMPO-CENTER  FOR  ADVANCED  STUDIES 
816  STATE  STREET 
PC  DRAWER  QQ 
SANTA  BARBARA  CA  93102 
B  GAMBILL 

D  AS  I  AC  (2l 

WARREN  S  KNAPP 

GEOPHYSICAL  INSTITUTE 
UNIVERSITY  OF  ALASKA 
FAIRBANKS.  AK  99  701 
TN  DAVIS 
NEAL  BROWN 
TECHNICAL  LABORATORY 

GTE  S  Y1.VANIA,  INC 
t-  cCTRONICS  SYSTEMS  GRF 
EASTERN  DIVISION 
77  A  STREET 
NEEDHAM. MA 

MARSHAL  CROSS 

II  T  RESEARCH  INSTITUTE 
10  WEST  35TH  STREET 
CHICAGO.  IL  60616 

TECHNICAL  LIBRARY 

UNIVERSITY  OF  ILLINOIS 

DEPARTMENT  OF  ELECTRICAL  ENGINEERING 
UP  BAN  A .  IL  61803 

AERONOMY  LABORATORY  !2> 

JOHNS  HOPKINS  UNIVERSITY 
APPLIED  PHYSICS  LABORATORY 
JOHNS  HOPKINS  ROAD 
t.  AuREl,  MD  20810 
J  NEWLAND 
PT  KOM1SKE 

-OCKHEED  MISSILES  &  SPACE  CO.  INC 
3251  HANOVER  STREET 
PALO  ALTO.  CA  94304 
EE  GAINES 
Wl  IMHOF  D, '52-12 
JB  REAGAN  0652-12 
RG  JOHNSON  D'52-12 

-..JWE-L  RESEARCH  FOUNDATION 
450  AIKEN  STREET 
lOVVELL.  MA  01854 

of;  bi  bl 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
LINCOLN  LABORATORY 
PO  BOX  73 

LEXINGTON,  MA  02173 
DAVE  WHITE 
jh  PANNELL  L  246 
DM  TOWLE 


MISS  ION  RESEARCH  CORPORATION 
735  STAI  E  STREET 
SANTA  BARBARA,  CA  93101 
R  HENDRICK 
F  FAJEN 
M  SCHEIBE 
J  GILBERT 
CL  LONGMIRE 

MITRE  CORPORATION 
PO  BOX  208 
BEDFORD,  MA  01730 
G  HARDING 

PACIFIC  SIERRA  RESEARCH  CORP 
1456  CLOVER  FIELD  BLVD 
SANTA  MONICA,  CA  90404 
EC  FIELD,  JR 

PENNSYLVANIA  STATE  UNIVERSITY 
IONOSPHERIC  RESEARCH  LABORATORY 
318  ELECTRICAL  ENGINEERING  EAST 
UNIVERSITY  PARK,  PA  16802 

IONOSPHERIC  RSCH  LAB  !2> 

R&D  ASSOCIATES 
PO  BOX  9695 

MARINA  DEL  REY,  CA  90291 
FORREST  GILMORE 
WILLIAM  J  KAR2AS 
PHYLLIS  GREpFINGER 
CARL  GREIFINGER 
A  ORY 

BRYAN  GABBARD 
RP  TURCO 
SAUL  ALTSCHULER 

RAND  CORPORATION 
1700  MAIN  STREET 
SANTA  MONICA,  CA  90406 

TECHNICAL  LIBRARY  (2) 

CULLEN  CRAIN 

SRI  INTERNATIONAL 
333  RAVENSWOOD  AVENUE 
MENLO  PARK,  CA  94025 
DONALD  NEILSON 
GEORGE  CARPENTER 
WG  CHETNUT 
JR  PETERSON 
GARY  PRICE 


STANFORD  UNIVERSITY 
RADIO  SCIENCE  LABORATORY 
STANFORD,  CA  94305 
RA  HELLIWELL 
FRASERSMITH 
J  KATSUFRAKlS 

TRW  DEFENSE  &  SPACE  SYS  CROUP 
ONE  SPACE  PARK 
REDONDO  BEACH.  CA  90. '78 
DIANA  DEr 

CALIFORNIA  INSTITUTE  Of  ft  WNOLO 1.  > 
JET  PROPULSION  I.ABOHATOR' 

4800  OAK  GROVF  DRIVE 
PASADENA.  CA  91  V  03 
ERNEST  K  SMITH 
(MAIL  CODE  ’44-B13' 
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